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Abstract 

The current measurements of the cosmic ray energy spectrum at ultra-high energies (E > 10 19 
eV) are characterized by large systematic errors and poor statistics. In addition, the experimental 



it is crucial that the extent to which these measurements disagree be well understood. Here we 



fN| . results of the two experiments with the largest published data sets, AGASA and HiRes, appear 

>. 

' to be inconsistent with each other, with AGASA seeing an unabated continuation of the energy 

spectrum even at energies beyond the GZK cutoff energy at 10 19 ' 6 eV. Given the importance of the 

\o 

related astrophysical questions regarding the unknown origin of these highly energetic particles, 

o 

QiV evaluate the consistency of the two measurements for the first time with a model-independent 

2 ■ method that accounts for the large statistical and systematic errors of current measurements. We 

> 

further compare the AGASA and HiRes spectra with the recently presented Auger spectrum. The 
method directly compares two measurements, bypassing the introduction of theoretical models 

X' 

5_i ■ for the shape of the energy spectrum. The inconsistency between the observations is expressed 

& . 

in terms of a Bayes Factor, a standard statistic defined as the ratio of a separate parent source 
hypothesis to a single parent source hypothesis. Application to the data shows that the two-parent 
hypothesis is disfavored. We expand the method to allow comparisons between an experimental 
flux and that predicted by any model. 
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I. INTRODUCTION 



By measuring the energy spectrum of cosmic rays above 10 19 eV we hope to shed some 
light on the yet unknown sites and acceleration mechanisms that can produce a particle 
of such energies. These measurements are intrinsically difficult, as the cosmic ray flux at 
ultra-high energies is low and we rely on earthbound detectors which cannot directly detect 
the primary particle. Rather, its properties are reconstructed by measuring the secondary 
particles of the extensive air shower produced when the primary enters Earth's atmosphere. 
Low statistics and large systematic errors in the energy determination are typical for these 
measurements, which should consequently be treated with care. 

The existing measurements of the energy spectrum above 10 19 eV have received much 
attention. A longstanding hypothesis predicts that the flux of the highest energy cosmic 
rays should be suppressed above 10 19,6 eV as cosmic rays from distant sources will interact 
with the cosmic microwave background via photopion production until their energy drops 
below this threshold energy Q, |J This so-called GZK suppression is in itself not a contro- 
versial prediction, but it has not been detected unambiguously, and several cosmic rays with 
substantially higher energies have been observed with various detectors over the years. 

The most recent results disfavoring the GZK hypothesis come from the Akeno Giant Air 
Shower Array (AGASA) collaboration, whose published energy spectrum shows no indica- 
tion of a high-energy suppression Q|. In contrast, the monocular-mode energy spectra 
measured by the High Resolution Fly's Eye detectors (HiRes 1 and HiRes 2) support the ex- 
istence of a GZK feature 0,0]. The current world data set is still small. The AGASA claim 
that the GZK supression is not observed is based on only 11 events above 10 20 eV. In the near 
future, the Pierre Auger Observatory, currently under construction in Malargue, Argentina, 
will dramatically increase the world data set. The Auger collaboration has published a first 
energy spectrum based on 1.5 years of data taken during construction ?]]. 

A second problem is that although typically not shown in plots of the energy spectrum, 
the errors on the energy determination are also large: 30 % at 3 x 10 19 eV and 25 % above 
10 20 eV in AGASA [4J, 30% in HiRes [J and 50% in Auger Q. Furthermore, a correct 
evaluation of the systematic errors is difficult. It should for example be noted that the 
systematic uncertainty for the AGASA energy scale includes a 10 % systematic uncertainty 
due to the hadronic model. This uncertainty is calculated by comparing different hadronic 
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event generators and denning the systematic uncertainty by their difference. This may not 
be a reliable estimate of this uncertainty, which could potentially be much larger and is, in 
any case, unknown at this point. Since the AGASA and HiRes experiments use radically 
different techniques to determine the primary energy - AGASA is a ground array and HiRes 
an air fluorescence detector - the errors have fundamentally different sources. 1 

The large statistical and systematic errors quoted by the experiments raise the question 
how much significance should be attached to the discrepancy between the HiRes and AGASA 
spectra. Several authors have addressed this question. In 0,13) the experimental results 
are compared individually to model predictions of the shape of the energy spectrum at GZK 
energies. Such an analysis requires assumptions on the nature of the cosmic ray sources, as 
the shape of the GZK suppression will of course depend on the source distribution. 

Before one asks if any of the current measurements can be used as evidence for or against 
a specific model for the origin of cosmic rays, one needs to answer the question of whether 
or not the two experiments actually disagree at all, given their uncertainties. This requires 
a method to compare the two spectra that is both model-independent and considers the 
uncertainties quoted by the experiments by accounting for the probability of any systematic 
shift in energy scale applied to the data. 

This paper presents a general spectrum comparison technique that for the first time ad- 
dresses these points. The technique naturally incorporates the Poisson errors in the observed 
fluxes, and can probabilistically treat systematic errors in the absolute energy scale. The 
technique is developed in Section 2 and applied to the energy spectra of AGASA jj], HiRes 
1 and HiRes 2 in monocular mode 0, Q and Auger [s| in Section 3. Section 4 summa- 
rizes the results and describes how the method can be expanded for testing all experimental 
measurements against a theoretical prediction. 



II. PROPOSED COMPARISON BETWEEN AGASA, HIRES AND AUGER 

The comparison between the AGASA and HiRes spectra can be quantified as the relative 
probability that the various fluxes came from a signal parent distribution versus two separate 

1 It should be noted that the published Auger spectrum has an energy scale that is calibrated with its 
fluorescence detector. If simulations are used to determine the shower energies from surface detector data 
alone, the energies are systematically higher by at least 250- 
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FIG. 1: The AGASA ||, Auger 



and HiRes monocular spectra 
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distributions. This probability ratio is quantified as a Bayes Factor (BF). Two parent fluxes 
can arise for any number of reasons, for example systematic mismeasurement by either 
experiment, or - however unlikely - a different primary flux at AGASA or HiRes. 

The calculation of the BF is a variation of the method given in [11] comparing the 
consistency of two Poisson parameters. The statistic is the ratio of the probability that 
the spectra arise from separate parent distributions over the probability of the alternative 
hypothesis that they are derived from a single parent distribution. The BF then indicates 
to what degree the separate parent hypothesis is favored. A widely accepted interpretation 
of the resulting value for the BF is shown in TablcQ] 

To perform the calculation, we divide each energy distribution into N bins. For the i th 
energy bin, the relative number of events is parametrized by the variable If the single- 
parent hypothesis is true, the relative exposures of the two experiments should determine 
the number of events expected from that single parent distribution. In the case of the two- 
parent hypothesis, a similar parameter, is invoked that can take any value for any bin 
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TABLE I: The interpretation of the separate parent flux hypothesis (SPFH) given the Bayes 
Factor (BF) Q- 



Bayes Factor (BF) 


Interpretation 


BF > 1 


SPFH supported 


10~s < BF < 1 


Minimal evidence against SPFH 


10 _1 < BF < 10~i 


Substantial evidence against SPFH 


1(T 2 < BF < io _1 


Strong evidence against SPFH 


BF < 1(T 2 


Decisive evidence against SPFH 



and is not constrained by the relative exposures of the experiments. The probabilities of 
both hypotheses are calculated assuming Poisson uncertainties for the fluxes in the i th bin, 
while marginalizing the energy scale uncertainties. 

We will derive the BF using the comparison of the AGASA and HiRes spectra as an 
example; however, in the end, we will calculate a separate BF for every combination of 
AGASA, HiRes 1, HiRes 2 and Auger spectra. The spectra for AGASA and HiRes at a 
given energy scale are defined as Aj and Hi, where i G [1,-ZV]. We define the percent shift 
in the AGASA and HiRes energy scales to be sa and sh, with the percent energy scale 



uncertainty to be a a and an, respectively. Following the derivation in [11 1. the true number 
of events measured by AGASA and HiRes in the i th energy bin is parameterized by 

Oi = fiVi (1) 

and 

hi = (l- fi)rh, (2) 

respectively, where rji is the total number of hypothesized events expected for both the 

AGASA and HiRes experiments. For the hypothesis where the two spectra come from a 

single parent distribution, 

j, (AGASA Exposure)i 

1 (AGASA Exposure)i + (HiRes Exposure)i 

In the case where the AGASA and HiRes experiments are measuring separate parent 

distributions, we parametrize the hypothesized events in a similar fashion. For the separate 



parent hypothesis, 




(4) 



and 



k = (i- m 



(5) 



where // is allowed to take any value and hence will be marginalized. 
Making use of Bayes' Theorem, the Bayes Factor is 



P(separate parent hypothesis\A, H) 




(6) 



JP(A,H\f, V )q(f, V )d V 



where, e.g., A = {Ai}. If one had a model with which to constrain A and H, it would 
be included in the prior, q(.). Any unknowns in the model could then be marginalized 
with r)i and j[ provided that the priors on r)i and // do not cancel. The denominator of 
Eq. (JHJ) assumes some values for /j calculated through Eq. (JHJ). The numerator marginalizes 
f- assuming any and all sets of By marginalizing over we are effectively asserting 
that the relative exposure quoted by the two experiments (defined by fa) does not hold. 
Therefore, some adjustment needs to be made to one or both of the exposures. Since we 
do not have any idea what that adjustment could be, we sum over every possible relative 
exposure (i.e. marginalize //). Note: if the difference in the spectra are due to some physical 
process, this would manifest itself in a change in the relative exposures. 
To introduce the energy scale uncertainties, we modify 



to account for the deviations of the true energy scale from the mean energy scale of AGASA 
and HiRes, respectively. A and H are made dependent on the energy scale, such that 
the values of Ai and Hi will depend on sa and s#, respectively. 1 If the numerator and 

1 Some may interpret this as changing the measurements, and Hi, based on the choice of energy scale. 
We can escape this interpretation of "changing measurements" if we imagine that, as the various terms for 
sa and s# are calculated, those events that fall inside the "energy window" between 10 196 and 10 21 eV 
are calculated according to Eq. ©. The probability distribution for those events whose true energies fall 
outside this window is flat (= 1). That is, we impose the arbitrary condition that if an event's true 



P(A, H\f, rj)q(f, n) H\f, r,, s A , s H ) 



(7) 
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denominator are properly evaluated in this form, they necessarily account for the probability 
of some observed fluxes, A and H, given the energy scale. As the deviations of the true 
energies from the measured energies are unknown, probability theory allows sa and s# to 
be marginalized such that Eq. (jUJ) becomes 

BF = 1 1 1 1 P(A H\f, 77, sa, s H )q(f, V, sa, s H )df'dr]ds A ds H 
1 1 1 P{A,H\f,r},SA,s H )q(f,r},SA,s H )drjds A ds H 

(8) 

The fluctuations in the fluxes are treated as Poisson, while the distributions of sa and 
sh are parametrized by Gaussians. Explicitly, the BF is 





[nf=i [/ / p nm {Ai)p (l - m mdfM ) 


\ N(s A ] 0, a A )N(sH] 0, a H )ds A ds H 


II\ 


[nf=i u Pf^Pii-f^mdvi]) 


\ N(s A ; 0, a A )N(s H ] 0, a H )ds A ds H 



(9) 



where P /i (M) denotes a Poisson distribution with mean /x and M measured events and 
N(x; fi, a) denotes a Gaussian with a deviation, x, from the mean, /x, and error a. By 
integrating 77^ within the product, we are effectively summing over every imaginable set of 
Poisson-distributed spectra. Each set of A and H, where A4 and Hi are dependent on the 
energy shift, is weighted by the probability of that shift. Since the size of the shift toward 
the true energy scale is unknown, these shifts are marginalized by integrating sa and s#. 
The numerator has an integration over //, effectively marginalizing over the cases where A 
and H come from different parent distributions. 

To get a sense of the behavior of the BF, we apply the method to a simple toy example. 
We generate two distributions, S and B, each with 20 bins; distribution S has 200 events, 
and distribution B has j3 x 200 events such that fi = where fi is identical for every bin. 
In effect, j3 defines the relative exposures of the two "experiments" S and B. The number 
of events in each bin of S and B fluctuates according to Poisson statistics. 

On top of these "statistical" fluctuations, we introduce a second fluctuation of distribution 
B, while the mean of the S distribution is kept flat. Each bin of the distribution B is evenly 
fluctuated within ±e, constraining the total events in B to 200 x (3, thereby making the 

energy is < 10 19 6 eV or > 10 21 eV, then nothing can be said about the predicted number of events in 
that particular bin. In this way, we can calculate the BF for a finite energy range, avoid the concept of a 
changing measurement based on a changing hypothesis and still shift the spectra relative to one another 
in a somewhat intuitive manner. 
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FIG. 2: BF as a function of e (the difference between the two parent distributions) and (3 where 
the BF is evaluated over two spectra with 200 and 200 x (3 events. 

distributions more different as e increases. Therefore, e quantifies the disagreement between 
S and B. Note that in contrast to the statistical fluctuations, these e fluctuations are not 
included in the BF and therefore deviations between S and B due to such fluctuations are 
'perceived' by the BF as differences in the spectra. We then plot the BF as a function of 
(5 and e in Fig.|21 

The BF distribution, given e and /, has substantial tails that cause large fluctuations in 
Fig. |21 even when averaged over hundreds of thousands of trials. Despite these fluctuations, 
we are still able to pick out trends in the BF as a function of (3 and e. For instance, the BF 
systematically increases with e (i.e. as the difference between the distributions increases). As 
(5 increases, the statistics for B increase. Therefore, it is more likely that a given difference 
between two distributions is due to a real difference in the parents rather than low statistics. 
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III. IMPLEMENTATION 



Due to practical limitations for the calculation of Eq. (|9"J). some choices must be made in 
the binning and shifting of the data. The calculation occurs over 14 bins of size logE/eV = 
0.1 from 10 19,6 eV to 10 210 eV. The energy bins are determined from existing data. The 
10 19,6 eV cut-off is motivated by a previous comparison in Q. The AGASA experiment has 
866 events above 10 19 eV and 72 above 10 19,6 eV. HiRes 1 and HiRes 2 have 403 and 95 events 
above 10 19 eV and 35 and 6 events above 10 19 ' 6 eV, respectively. Auger has 444 events above 
10 19 eV and 17 above 10 19 - 6 eV. 

After integrating over // and rji, Eq. (jUJ) reduces to 



BF 



II 


\u N 1 

lli=l Ai+Hi+1 


N(s A ; 0, a A )N(s H ; 0, a H )ds A ds H 


II 


lli=lJi \ L JiJ AilHil 


N(s A ; 0, <r A )N(s H ; 0, a H )ds A ds H 



(10) 



Since one cannot split a spectrum into an infinite number of bins, one must make finite 
shifts in the two energy spectra. The chosen shifts are made ±5cr from the mean energy, 
and are incremented in log-E/eV = 0.1 steps, weighted by the Gaussians integrated over 
the energy bin. As various energy scales sa and sh, are evaluated, the fluxes in the 14 bins 
change. The shifts s A and sh are evaluated down (up) to logE/eV = —0.5 (logE/eV = 0.5). 
For instance, energy bins below 10 19 6 eV shift in and out of the set of 14 bins depending on 
the values of sa and sh- In the case of the HiRes 2 spectrum where the data are quoted 
in increments of logE/eV = 0.2, the events are divided in proportion to the absolute size 
of the logE/eV = 0.1 bin. Meanwhile, the relative exposures of AGASA and HiRes (/) 
remain unchanged. 

Fig. 01 shows the BF as a function of the percent energy uncertainties for all combinations 
of AGASA, HiRes 1, HiRes 2 and Auger spectra. The distributions are non-trivial to 
understand as changing the energy scale does not amount to a simple shift in the number of 
events from one energy bin to another. Rather, in HiRes, the shape of the spectra change 
as the energy scale changes because the exposure remains fixed as a function of energy. 
Further, the probability that the spectra agree for a given energy scale is convoluted with 
the probability that that energy scale is the true energy scale for the detector. Meanwhile, 
the lower energy bins fall in and out of the calculation as the energy scales change in the 
course of the calculation. 
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If we compare AGASA and HiRes with 30% energy scale uncertainties for both experi- 
ments, the BF = 0.71 and BF = 0.04 for HiRes 1 and HiRes 2, respectively. According to 
TableEJ the BF for the HiRes 1 comparison corresponds to minimal evidence against the 
hypothesis that AGASA and HiRes are derived from separate distributions. The latter value 
for HiRes 2 corresponds to strong evidence against the separate source hypothesis. If we 
compare the Auger spectrum with 50 % energy scale uncertainties with the HiRes spectra 
with 30% energy scale uncertainties, the BF = 0.54 and BF = 0.85 for HiRes 1 and HiRes 
2, respectively, corresponding to minimal evidence against the separate source hypotheses 
(i.e. no support for the two-source hypothesis). Finally, if we compare the Auger spectrum 
with 50 % energy scale uncertainties with the AGASA spectrum with 30% energy scale un- 
certainties we obtain a BF = 0.74 again corresponding to minimal evidence against the 
separate source hypothesis. 

In order to verify that the BF we obtain from the experimental data is indeed consistent 
with the single source hypothesis, we compare the actual BF for the spectra from experiment 
X and Y to simulations of the BF generated under the assumption that they are generated 
from a single parent. 

To calculate a simulated BF for the spectra from experiments X and Y, we assume 
that the true spectrum is a spectrum that is an average, weighted by the exposure, of the 
spectra measured from experiments X and Y. This average spectrum is then scaled to the 
corresponding exposures for X and Y. We then randomly fluctuate both spectra according 
to Poisson statistics and calculate the BF. This procedure is repeated many times. The 
resulting distribution of the BF for two spectra generated assuming a single parent is then 
compared to the experimental value of the BF for the two measured spectra. Fig.H] shows 
the results for a comparison of all three experiments with each other. 

Comparing the BF's measured for data (denoted by a red line) and those calculated in 
Fig.Hl we see that, in general, there is a substantial probability of obtaining a BF that 
is larger than the BF's calculated for the data. That is, a substantial number of samples 
that are derived from a single source look more like they were derived from two sources 
than the real AGASA, HiRes and Auger comparisons. The value of the BF tells us that 
the hypothesis that the AGASA and HiRes energy spectra are consistent with a single 
distribution is favored over the hypothesis that they are derived from separate distributions. 
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FIG. 3: The Bayes Factor (BF) as a function of the percent energy scale uncertainties for the 
compared spectra. 

IV. DISCUSSION 

This result is similar to those obtained for the AGASA and HiRes spectra in (t| 10] by 
a rather different method. Our method directly compares two experimental results rather 
than comparing each to a theoretical prediction. It also treats systematic errors correctly 
by weighting possible systematic shifts in the energy distribution by the probability of that 
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FIG. 4: The distribution of the BF calculated for an averaged spectrum from two experiments. 
The red lines denote the BF's measured for data. 



shift. 

This result does not say anything about what the specific (in) consistencies would be of ei- 
ther experiment given a theory. However, our results indicate that there is room for a theory 
that agrees with each pair of spectra. As the statistical power of the world data improves, 
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the methods discussed here will provide a convenient framework to compare experimental 
results with one another and with theoretical models. 

The Auger experiment is currently the only operating ultrahigh energy cosmic ray detec- 
tor. The Auger data set will increase dramatically over the next years. Figure|5] shows the 
projected value of the BF between Auger and the other experiments as a function of the 
fraction of the current Auger exposure. The BF is calculated assuming the current spectra, 
the existing 30 % energy uncertainties in AGASA and HiRes and the projected 15 % energy 
uncertainties in the Auger spectrum. In the limit where the exposure for Auger becomes 
large, the BF decreases (increasing support for the one-source hypothesis). That is, the 
statistical and energy scale uncertainties in HiRes and AGASA are so large that, barring a 
significant change in the Auger spectrum, one will never be able to tell whether or not the 
Auger and AGASA (HiRes) spectra are derived from different parent distributions provided 
that we limit our analysis to events above 10 19 ' 6 eV. Note, however, that if one fits a model 
to one of the experimental spectra, our result does not imply consistency of the model with 
the other experiment. 2 

This project is supported by the National Science Foundation under contract NSF-PHY- 



2 In order to compare a model with both experiments, one must compare the model to both data sets 
simultaneously. The above method can be modified to provide a statistic analogous to a x 2 probability 
that could measure the consistency of theory given all experimental results. In this case, Bayes' theorem 
is used to obtain 

_ P(A,H\i) q (i) 
P{t\A,H) - ^ theo „ es P(A,H\t)q(t)' (U) 

where, because of the definition of / (/'), each ti is the total number of events expected in the i th bin for 
both experiments, and, likewise, each U is the number of events expected from any physically meaningful 
theory. Using the experimental Poisson and Gaussian errors, 



P(t\A,H) = 



II 


nf =1 




N(s A ;0, a A )N(s H ; 0, a H )ds A ds H 


Ai'.Hil 




II 


\v\ N 1 

lli=l Ai+Ht + l 


N(s A ; 0, <r A )N(s H ; 0, a H )ds A ds H 



(12) 



where all theories are effectively marginalized by integrating within the product. To obtain a x 2 -like 
probability, one would calculate P(i\A,H), and then calculate many P(t'\A, H)'s where each if. is the 
Poisson fluctuated tj. The fraction of events where P^'^A^H) < P(ti\A,H) would be the degree of 
confidence in the theory. 
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FIG. 5: The projected BF between the current Auger spectrum and the AGASA and HiRes 
spectra as a function of the fraction of the current Auger exposure. The calculations assume 30 % 
energy uncertainties in AGASA and HiRes and the projected 15 % energy uncertainties in the 
Auger spectrum. 
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